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ABSTRACT 

We present 2. 5- dimensional time-dependent simulations of nonrelativistic and 
nonradiative outflows from sinusoidally perturbed Keplerian accretion disks. A 
sinusoidal perturbation is introduced in the velocity of the gas ejected from the 
surface of the disk into a cold corona. In the simulations, the disk is a flxed 
boundary from which the gas is ejected with a pulsed velocity. The maximum 
value of this velocity is taken to be a thousandth of the local Keplerian disk 
velocity. It was found that for large periods, the structures in the jet tend to 
fragment into smaller substructures. For small values of the period, the structures 
tend to dissipate, while for medium values of the period, they tend to persist. 

Subject headings: ISM: jets and outflows — galaxies: jets — accretion, accretion 
disks— MHD. 

1. Introduction 

In general, the production of astrophysical jets has been primarily explained by two 
possible mechanisms: magnetic forces associated with an accretion disk (Konigl & Pudritz 
1999) and the interaction between an accretion disk and a magnetized star, in the case of 
protostellar jets (Shu et al. 1999). The flrst mechanism is a disk- wind driven model, where 
the gas is centrifugally flung out from the surface of a Keplerian disk when field lines thread 
the disk at an angle of 60° or less with respect to the surface. When the gas reaches the 
local Alfven velocity, the ( Jp x B^) term becomes dominant, providing the coUimation for 
the outflow. The second mechanism is the X-wind model, where the wind emanates near 
the region where the disk corotates with the central magnetic star. It is assumed that 
some magnetospheric lines, connecting the central star and the disk, are partially open. 



-3- 



Excess angular momentum brought in by the disk is extracted by the magnetocentrifugally 
driven wind. The wind emanates from the inner edge of the disk and flows out along 
the open magnetic field lines. Since the first mechanism seeks to explain astrophysical jets 
independently of the detailed nature of the central object, from active galactic nuclei (AGNs) 
to young stellar objects (YSOs), it has the possibihty of being a universal model. 

Recent theories and models of astrophysical jets have mainly tried to provide a natural 
explanation for common features. These include the episodic outbursts and knots that are 
observed. 

Numerical simulations of axisymmetric magnetized flows from accrection disks have 
previously been made using different models. Shibata & Uchida (1986) and Stone & Norman 
(1994) followed the internal dynamics of the disk. Bell & Lucek (1995) and Ustyugova et 
al. (1995) treat the disk as a boundary condition at the base of the wind, without following 
its internal dynamics, and assume it to be in Keplerian equilibrium. Following the latter 
model, Ouyed & Pudritz (1997a, b; 1999) made 2.5-dimensional, high-resolution, numerical 
magnetohydrodynamic (MHD) simulations of the onset and coUimation of outflows from the 
surface of a Keplerian accretion disk around a central object. Using two different topologies 
for the magnetic field that threads the disk and the corona, they obtained steady or episodic 
jets, depending on the values of the constant (in time) ejection velocity from the disk to 
the corona. For the smallest values of the ejection velocity, their results showed episodic 
production of knots near the surface of the disk which propagate along the jet axis. They 
argued that these knots arose due to MHD shocks. 

Krasnopolsky et al. (1999) made time-dependent, axisymmetric simulations of the mag- 
neto centrifugal model jet formation. Their system consisted of a Keplerian disk and a central 
object. The main subtlety resides in the treatment of the disk- wind boundary, the disk sur- 
face. Some, but not all, of the flow variables can be fixed on this boundary. The boundary 
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conditions imposed on the disk surface were limited to those necessary and sufficient to take 
into account information propagating upstream from the fast and Alfven critical surfaces, 
avoiding overdetermination of the flow and the production of impulsive accelerations. They 
obtained a cold and steady jet, launched smoothly from the Keplerian disk and coUimated 
by the magnetocentrifugal mechanism. The coUimation is observed both in the shape of the 
field lines and in the density profiles, which become cylindrical, and thus jetlike, at large 
distances along the rotation axis. The structure of the flow is insensitive to the density or 
flow speed at the injection surface, as long as the mass flux is kept constant. As explained 
by the authors, steady-state jets obtained in earlier studies, using an essentially incomplete 
treatment of the disk boundary conditions, are qualitatively unchanged when this deficiency 
is rectified. It remains to be shown, however, whether the same is true for nonsteady jets. 

Frank et al. (2000) made a series of simulations intended to address the issue of MHD jet 
propagation. Whereas some studies have used ad hoc initial conditions, they injected flows 
into a computational grid derived from models of coUimated jets driven by magnetocentrifu- 
gal launching. The initial configurations in the jet were taken directly from the solution of 
the force balance perpendicular (the Grad-Shafranov equation) and parallel (the Bernoulli 
equation) to the magnetic surfaces generated by a magnetized rotator. Their simulations 
followed the evolution of jets composed of helical fields, embedded in hypersonic plasmas 
whose density and velocity varied with radius. They studied the propagation of jets driven 
by: (1) a purely Keplerian rotator; and (2) a Keplerian rotator with a sub-Keplerian bound- 
ary layer. Both rotators are exterior to a solid body rotator. Studying the adiabatic case for 
the last model (their multicomponent model) they obtained a jet in which strong instabilities 
developed in the core of the jet. Once the core is exposed, periodic pinches can be seen in 
the beam. Their stability analysis indicates that the origin of these periodic pinches could 
be current driven instabilities. The pinching instabilities on the axis have a wavelength of 
approximately half of the jet radius. However, the beam also has a second characteristic 
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wavelength which runs along the surface of the core. This feature, which appears as an 
envelope encompassing the shorter wavelength modes, has a wavelength of approximately 
three jet radii. As explained by the authors, these results suggest that these instabilities 
could be at the origin of the knotty structure of a large number of jets as seen, for example, 
in HL Tau, HHl, HH30, and HH34. 

Time-dependent MHD calculations of a modified version of the X-wind model, consid- 
ering the star-disk field interaction, shows twisting of closed stellar field hues, current sheet 
formation, reconnection, and ejection of magnetic island with the plasma heated to < 10^ 
K, consistent with X-ray flare observations (Hayashi et al. 1996; Goodson et al. 1997). 

All the above previous studies assumed constant velocity injection from the disk. Here, 
we extend the previous investigations and study the effect of non-steady injection, in partic- 
ular, periodic injection. 

Adopting the disk as a boundary condition, we investigated how periodic perturbations 
at the disk surface can alter the formation and the propagation of the magneto-centrifugal 
jet. In the FU-Orionis phase of a protostar, for example, we know that there is a periodic 
ejection of material observed in the X-ray range. We do not include the physics of the 
disk that could generate a pulsed gas ejection. We assume that there is a periodic mass 
loading at the base of the corona of the disks, which introduces a periodic ejection of gas 
from the disk surface. Whereas Ouyed & Pudritz (1997a, b; 1999) used a constant ejection 
velocity from the disk, we study the temporal dependence of this velocity. We made 2.5- 
dimensional numerical simulations using the ZEUS-3D code (Stone & Norman 1992a,b) in 
the axisymmetry option. A sinusoidal perturbation was introduced in the ejection velocity 
as a temporal function. Various values for the periods in the range from 10 to 80 rotation 
periods of the innermost radius of the accretion disk (r^), were used in the simulations. The 
maximum value for the perturbation was chosen. This value is roughly equal to that which 
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would produce a steady state jet if the gas were continuously ejected from the disk (Ouyed 
& Pudritz 1997a, hereafter OP97a). The luminosity of the protostellar T Tauri stars (which 
have accrection disks) fluctuates rapidly with large amplitude. Similarly, quasars, which are 
also assumed to have accretion disks, also have luminosities which fluctuate rapidly with 
large amplitude. Therefore, we introduced in our simulations velocity perturbations that 
are a large percentage of the injection velocity. A region near the central object, SOrj in the 
2;-direction and 20rj in the r-direction, was simulated. 

Using different values for the period of the sinusoidal perturbation, we obtained the 
formation of regularly spaced structures along the jet axis. For small values of the period, 
the structures tend to dissipate along the jet axis; for medium values, they tend to persist 
and for large values, tend to fragment into smaller substructures. 

In §2 the basic equations that are solved are presented. The model is described in §3 
and the numerical results are presented and discussed in §4. Finally, in §5 we present a 
summary and the final conclusions. 



2. Basic equations 

Adopting cylindrical coordinates (r, 4>, z), we placed the central object at the origin 
and took the z-axis to be perpendicular to the disk. The surface of the disk thus lies in the 
z — plane of our coordinate system. 

The equations that describe the conservation of mass, momentum, energy and magnetic 
flux are: 

^ + V-{pv) = 0, (1) 
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dv 
'dt 



{v-V)v 



+ V{p + pa)+ pV$ - i X B = 0, 



(2) 



de 
di 



+ p{V ■v) = 0, 



(3) 



dB 

'dt 



- V X (v X B) = 0, 



(4) 



V • B = 0, (5) 

with p the gas density, v the velocity, B the magnetic field, p the gas pressure, e the 
internal energy, j — V x B/Att the electric current density, and $ = —GM/ (r^ + 2;^)^^^ the 
gravitational potential. 

A polytropic gas p — Kp^ with the polytropic index 7 = 5/3 is assumed. Thus, ignoring 
radiative transfer effects, we do not solve the energy equation (3). We also assume that pa is 
an Alfvenic turbulent pressure. This provides, together with the thermal pressure, the main 
support against gravity for the cold corona. Here pa — p/ Pti where (5t is the ratio of the 
gas pressure to the Alfvenic turbulent pressure, assumed to be constant (OP97a). 

All parameters are given in units of their values at the innermost radius of the accretion 
disk Ti. Thus, the normalized equation of motion (eq. [2]) is 



f + (.'.v>' 



5i p' Si/3i p' 

where the dimensionless variables are: r' — r/ri, z' — z/ri, v' — v/vK,i, p' — p/pi, p' — pjpii 

I 1 /2 

B' = B / Bi, $' = —1/ (r'^ + z"^^ and V = rjV. The parameter 5i = Pivj^^/pi is the ratio 
of the Keplerian to the thermal energy density and /3j = Snpi/B^ is the ratio of the gas to 
the magnetic pressure in the corona at rj. At rj, the Keplerian speed is VK,i — (GM/ri)^^'^ 
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and Bi is the poloidal field at rj. For convenience, we drop the primes and refer only to the 
normalized variables. Time is given in units of tj = ri/vK,i, so that the dimensionless time is 

r - p (7) 



3. The Model 

In this model, the disk is a fixed boundary condition in pressure equilibrium with the 
corona. At the center there is a point mass that represents a star or a black hole. An initial 
potential (zero current) poloidal field threads the disk and the corona. 

The computational domain defines the spatial region where the MHD equations are time 
evolved. ZEUS-3D, a grid code, divides the computational domain into cells, called "active" 
zones. Finite differencing the evolution near the boundaries of the grid requires that values 
for the dependent variables be specified beyond the computational domain. Therefore, at 
each boundary, "ghost" zones arc added. Values for the dependent variables in the ghost 
zones are specified, using boundary conditions appropriate to the geometry and the physics 
of the problem being solved. Thus, the evolution equations are not solved for the ghost 
zones. 

We use inflow boundary conditions at the disk surface and outflow boundary conditions 
on the remaining boundaries. Reflecting boundaries are used along the axis of symmetry, 
which coincides with the disk axis. The normal component (radial component) of the velocity 
and the magnetic field are set to zero on the boundary and reflected into the "ghost" zones. 
In addition, the 3-component (toroidal component) of the velocity and magnetic field is also 
inverted along the axis of symmetry into the "ghost" zones. The initial setup and boundaries 
are shown in Figure 1 (Fig. 2 of OP97a). 



EDITOR: PLACE FIGURE 1 HERE. 
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In the corona, the initial toroidal field is zero. Magnetostatic equilibrium is numerically 
maintained to assure that an outflow can only arise due to the magnetocentrifugal mech- 
anism. Resolving equation (6) for magnetostatic equilibrium, we obtain the expression for 
the coronal density distribution, 

p=(\j^?^^^''' ^ (8) 

which is used in our numerical simulations. 

Since the surface of the accretion disk is a fixed boundary condition, there is no back 
reaction by the jet. We are, thus, time evaluating the coronal region (2; > 0) only where 
the jet arises. However, we fix the density at the base of the corona (first column of cells), 
assuming that the disk is a suitable dynamic source of matter. The density distribution of 
the disk surface has the form 

Pd = mr-^'\ (9) 
where 77^ = Pd/ Pq is the ratio between the disk density (p^) and the density at the base of 
the corona (po = p{z = 0)). 

The corona and the disk are threaded by the poloidal magnetic field in a potential 
configuration (J = 0), described by the stream function, 

= \lr'' + {zd + zf - {Zd + z) , (10) 



where z^, is the dimensionless disk thickness. The field components are given as 

= — and = (11) 

r or r oz 

Due to the torsion of the poloidal field at the disk surface, a toroidal magnetic field is 
produced at the base of the corona, which could propagate into the disk. We thus expect 
that the disk has a nonzero toroidal field, 

B, = ^, (12) 
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where the input parameter /li — B^i/Bi is the ratio of the toroidal to poloidal magnetic field 
at Tj . There are three important initial timescales defined at the disk's surface: the Kepler 
time tk, the Alfven time ta and the magnetic braking time tb- Thus, we have (e.g., eq. 
[4.64] of OP97a) 

[Tk, Ta, Tb) = (i, y/5Aj2, mVW^) ■ (13) 

To estimate the turbulent pressure needed to support the corona for a given gas pressure, 
we use the equation for the ratio of the turbulent to magnetic pressure /3t 

/3t = 7/((7-1)^^-7), (14) 

where, as in equation (6), 5i is the ratio of the Kepler to thermal energy density. The gas is 
injected from the disk into the corona at a velocity Vp — fyVxBp/ Bp, where vk — {^/f)^^'^ 
is the local Keplerian velocity at the disk surface (for r > 1)^, Bp is the local poloidal 
magnetic field that threads the disk, and is the temporal function that describes the 
temporal dependence of the injection velocity. In our simulation we introduced a pulsed 
injection velocity with a sinusoidal form, 

fv = Vinjo [1 + sin {2'kt/T)] , (15) 

where T is the period of the perturbation, vinjo is an input parameter and r is given by 
equation (7). 

We, thus, have six free parameters in our simulation: (3i, rji, /li, Vinjo and T. 

^For the radial region where the disk does not exist (r < 1 (i.e., r < n)) we set Vp to 
zero. 
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4. Numerical results and discussion 

We used the following parameters to investigate the time evolution of a jet from a 
Keplerian accretion disk with pulsed gas injection (similarly to OP97a): 

{5i, f3i, rji, i^i, Vinjo) = (100.0, 1.0, 100.0, 1.0, 10-^) , (16) 

implying the initial timescales 

{tk, ta, tb) = (1.0, 7.07, 707.1) . (17) 

Five simulations were made. In simulation A, there was no perturbation; in simulations B 
to E, we introduced perturbations with various values for the period T. The simulations 
were made in the domain {z,r) = (80,20) with a resolution of (120,45) for 400 time units 
(< Tb)- Along the z-axis, 30 zones span a distance of lOr^ outwards from the disk surface. 
An additional 90 ratioed zones span an additional distance of 70ri, with each subsequent 
zone increasing in size by a factor of 1.017. Twenty uniform zones span an initial five disk 
radius in the radial direction. Outside of these uniform grid zones, 25 zones are ratioed, with 
each subsequent zone increasing in size by a factor 1.065. 

Figure 2 shows the initial conditions in the atmosphere. The density distribution (eq. 
[8]) is shown at the top and the initial magnetic (potential) configuration, at the bottom of 
the figure. Initially, the velocity and the toroidal magnetic field are zero. 

EDITOR: PLACE FIGURE 2 HERE. 

The results of the simulations are shown in a series of figures that are described below. 
Figures 3 to 7 show the distributions of density, poloidal velocity (vectors), poloidal and 
toroidal magnetic fields and toroidal velocity at 400 time units for the simulations A to E, 
respectively. The contour lines are isocontours in the r — z plane. 
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In simulation A (Fig. 3), no perturbation was introduced in the gas injection. Results 
in good agreement with those of OP97a were obtained. Prom the density distribution, shown 
at the top of this figure, we observe the formation of a highly coUimated steady state jet 
along the z-axis. The coUimation is due to the strong radial pinch force of a dominant 
toroidal field (fourth panel from top to bottom) that propagates into the corona. We obtain 
a maximum value for the poloidal velocity (second panel from top to bottom) on the order 
of the Kepler ian velocity at r^. 

EDITOR: PLACE PIGURE 3 HERE. 

A sinusoidal perturbation in the gas injection velocity (eq. [15]) at the disk surface was 
introduced in the other four simulations. As we are interested in studying the behaviour of 

the jet subject to this perturbation at its base, we chose a progressive range for the period 
T. In simulations B to E, we use T = 10, 20, 40, 80, respectively. 

Along the jet axis in simulation B for T = 10 (Fig. 4), we observe the formation of 
regularly spaced structures up to about z = 50rj. The distance of separation between two 
consecutive structures is about Sr^ (i.e., ~ {T/2) r^). The axial extent of the regularly spaced 
structures advance with the flow speed until it reaches the limiting distance of about 50rj. 
This pattern persists even for times greater than 400 time units. 

EDITOR: PLACE PIGURE 4 HERE. 

Our results were obtained using ratioed zones that allow for a smaller number of cal- 
culated zones and, therefore, a poorer resolution as compared with OP97a. To verify the 
sensitivity of our results to this poorer resolution, we made simulation C using only uniform 
zones at a resolution close to that used in the uniform zone region in the poorer resolution 
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case. This simulation was run in the domain {z, r) — (80, 20) with the resolution of (250, 80) 
uniform zones. For the others cases (D - E), we used the poorer resolution. 

Regularly spaced structures along the whole jet (Fig. 5) are seen in simulation C, where 
T = 20. These structures persist during the entire simulation, even for times greater than 
400 units. The distance of separation between two consecutive structures is about llrj, or 
roughly (T/2) We obtained these same results when using ratioed zones in the poorer 
resolution case. 

EDITOR: PLACE FIGURE 5 HERE. 

We observe the formation of regularly spaced structures along the entire jet as well 
as their fragmentation in simulation D for T = 40. Few minor structures are observed 
between the major structures (Fig. 6: top). Fragmentation of the structures can be seen in 
the isocontours of the toroidal velocity in the bottom panel of Figure 6. These structures 
persist during the entire simulation, even for times greater than 400 units. The distance of 
separation between two major structures is about 22rj, i.e., roughly (T/2) r^. 

EDITOR: PLACE FIGURE 6 HERE. 

In simulation E for T = 80, we observe the formation of regularly spaced structures 
along the entire jet (Fig. 7: top), with more pronounced fragmentation clearly seen in the 
plot of the toroidal velocity (Fig. 7: bottom). These structures persist during the entire 
simulation, even for times greater than 400 units. The distance of separation between two 
major structures is about 40ri, i.e., ~ (T/2)ri. 



EDITOR: PLACE FIGURE 7 HERE. 
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The axial distributions, for all the simulations, of density, axial velocity, toroidal velocity 
and toroidal magnetic field along a cut at r = Sr^, parallel to the disk's axis at 400 time 
units, are shown in Figure (Fig. 8). These quantities are normalized to their maximum 
values along this cut. We note from the axial velocity distribution (dotted hue), that the 
flow is rapidly accelerated along the field hne from the disk to roughly z — 15ri. After this 
point there is a "Hubble flow" constant acceleration, with Vz proportional to z for simulation 
A. However, for simulations B to E this constant acceleration character is modified, reflecting 
the pulsed nature of the jet. From the toroidal magnetic field distribution (dashed-dotted 
line) we observe the effect of torsion of the poloidal field lines due to the rotation of the disk. 
The regular structures, in simulations B - E, have their minimum mass concentration (solid 
line) at the maximum values of toroidal magnetic field (dashed-dotted line), refiecting jet 
constriction by this field. The shift between the toroidal magnetic field and the density in 
the pinched structures is a general result of long standing (e.g., Chan & Henriksen (1980)) 
and it has been associated with the periodic oscillations caused by the pinch instability due 
to the toroidal magnetic field. 



For simulations A to E at r = 400, the axial and radial mass fluxes are shown in Figure 
9. The mass flux is deflned as 



where A is the area. The axial mass flux^ refers to the quantity of mass per unit time 
which passes through a circular area A — irr"^ and is calculated for each z, using r — 20rj; 

^The computation of the axial mass flux involves computing the product of Vz and the 
area of each annulus for each radial zone [7r(r|^]^ — r|)], where j is an active zone and rj refers 



EDITOR: PLACE FIGURE 8 HERE. 




(18) 
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the radial mass flux^ refers to the quantity of mass per unit time which passes through the 
area A = 2Trrz and is shown for each z, using r = 20rj. These fluxes, given in units of 
TJii = 2TTrlpiVKi-i are important for the detection of the decoUimation of the jet. Their values 
agree with the results of OP97a. At the time r = 400, the axial fluxes at 2; = 80^ are 
comparable to those at 2; ~ 0; the small difference owing to decoUimation of the jet. As the 
axial mass flux decreases, the radial flux can be seen to increase by roughly the same amount. 
We can see the pulsed nature of the jet (simulations B to E) from the axial mass fluxes (solid 
lines) which oscillate around their average values. The radial mass fluxes in simulations B 
to E are similar to those of simulation A (without perturbation). In simulation B, the axial 
mass flux curve is not very modified by the density structures, as compared with the case 
for which there is no perturbation (simulation A); in simulations C to D, the axial mass flux 
curves are strongly modified by the density structures. Their average values, however, do 
not change. 

EDITOR: PLACE FIGURE 9 HERE. 

Figure 10 shows the sonic, Alfvenic and magnetosonic Mach numbers for all the sim- 
ulations at T = 400, in a cut at r = hvi along the z-axis. In the computation of the sonic 
Mach numbers, we use only the thermal pressure p in the sound speed. To compute the 
magnetosonic Mach number, we use the effective sound speed which includes the thermal 

to the position of the face or edge of the zone. At each the axial mass flux is computed 
from the sum of the axial mass flux contributions of all 45 radial zones. 

^The radial mass flux at each z involves computing the product of Vj. and the area of the 
annulus for a given radial zone \2T:r j{zk+i — Zk)], where j and k deflne each active zone, and Vj 
and Zk refers to the face/edge position of the zone, along the r and z direction, respectively. 
In particular, we computed the radial flux at each z using r = 20ri. 



-16- 



and the Alfvenic turbulent pressure, pa^- We note in Figure 10 that the jet is already su- 
personic and super- Alfvenic near the base of the corona and super-magnetosonic at z ~ 4ri. 
Only where density structures exist are the Mach numbers (Fig. 10) modified. Their median 
values, however, remain roughly the same. The maximum values of Alfvenic Mach num- 
bers (dashed-dotted line) are due to the decreased axial magnetic fields. Negative values 
correspond to a reversal of the axial magnetic field. 

EDITOR: PLACE FIGURE 10 HERE. 

5. Summary and Conclusions 

We made numerical 2.5-dimensional simulations of astrophysical jets produced from 
Keplerian disks subject to a sinusoidal perturbation of the gas injection velocity from the 
disk into the corona. Ratioed zones were used in order to obtain better resolution near the 
region of the generation of the jet with a poorer resolution outside. 

Five simulations were made with different values of the period of the sinusoidal pertur- 
bation. No perturbation was introduced in simulation A. In simulations B to E, we used 
periods T — 10, 20, 40 and 80, respectively, in units of ti — Ti/vxi, where n is the innermost 
radius of the disk and Vxi, its Keplerian velocity. 

For simulation A, where there was no perturbation, we obtained a steady state highly 
collimated jet with a maximum axial velocity of the order of Vxi after a time r = 400. In 
a cut at r = 5rj parallel to the disk's axis, the jet is supersonic and already super- Alfvenic 

^The magnetosonic speed is defined as (c^ +v\^)^^'^ where Cg — (7(p + Pa) / pY^'^ is the 
effective sound speed and vaz — B/ (Anp)^^^ is the axial Alfven velocity. 
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near the base of the corona and super- magnetosonic at 2; ~ 4ri. These results agree with 
those of OP97a, who used better resolution. 

Introducing a sinusoidal perturbation in the gas injection velocity, the results common, 
to simulations B to E, are: 

1. The formation of regularly spaced structures with a separation distance of roughly 
(T/2) Tj along the jet axis. 

2. The toroidal magnetic field is out of phase (~ 90°) with the density curve along the 
jet axis. 

3. The axial mass fluxes are modified by the periodic density structures, although the 
average remains roughly the same, as compared with the same curve for the case of no 
perturbation (simulation A). The axial fluxes at z = SOrj are comparable to the axial fluxes 
at 2; ~ 0. 

4. The radial mass fluxes are similar to those of simulation A (where there was no 
perturbation). However, differences arise when we examine a cut at r = Sr^. There is 
oscillatory behaviour in some of the simulations at this cut. 

5. In a cut at r = 5rj the Mach number distributions are modified only where density 
structures exist. The median values, however, remain roughly the same, and are similar to 
simulation A. 

The results that are different in some of these simulations are: 

a. In simulation B {T — 10) the structures persist only up to an axial distance of 
z ~ SOrj. 

b. In simulations D {T — 40) and E {T — 80), fragmentation of the structures occur. 
These simulations, introducing periodic perturbations in the accretion disk, produced 
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coUimated jets with a wavelength of {T/2)ri, where T is the temporal period of the per- 
turbation. The creation of structures with a specific distance separation was determined 
by the period T of the gas ejection from the disk surface. For example, we obtained a jet 
with a distance separation of about llr^ associated with a period of 20 time units for the 
perturbation. However, for different values of T, an interesting effect occurs: the structures 
tend to dissipate at large axial distances for smaller values of T and to fragment into smaller 
structures for greater values of T. 

Using this same numerical approach, Ouycd & Pudritz(1999) obtained episodic jets, 
using steady gas injection from the disk with velocities which were smaller than lO"^^;^^^ 
which we used. For injection velocities equal to this value, they obtained a highly coUimated 
steady state jet^ with features similar to ours, with no perturbation in the velocity field 
of the disk. When we introduced perturbations with specific features in our steady state 
jet, our results showed the tendency to produce regular structures with separation distances 
~ llrj, similar to the case where OP99 used a smaller velocity of 10~^VK,i- 

Gardiner & Frank (2000), using a magnetocentrifugal approach for the jet launching, 
studied an analytical one-dimensional model, in which they investigated the behaviour of 
a jet with a variable output speed. They inferred that a pulsing jet with an embedded 
helical magnetic field will inevitably develop a periodic field structure, consisting of rarefied 
regions (where the helix is "combed out" , resulting in a more poloidal geometry) alternating 
with denser, toroidally dominated knots. To confirm these predictions, they performed 
axisymmetric (2. 5- dimensional) numerical simulations of a radiative MHD jet. They used 
conditions appropriate for protostellar jets and introduced a pulsed velocity with a sinusoidal 

^They obtained steady state jets with velocities equal to or greater than 5 x 10~^VK,i, as 
well as episodic jets with different distances of separation for smaller velocities, using the 
same free parameters that we used. 
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form at the jet nozzle. Their results showed that pulsing jets with initially helical fields will 
evolve into an alternating poloidal-toroidal-dominated field geometry. The strongest toroidal 
field are confined to the dense knots, while predominantly poloidal fields are in the rarefied 
regions. 

Both our results and those of Gardiner & Frank (2000) show that episodic jets are 
produced by introducing a pulsed mass load at the base of the jet. In particular, in our 
simulations, as well as in the study made by Ouyed & Pudritz (1997b), knot formation is 
associated with an anticorrelation of density and toroidal field strength. However, it should 
be noted that our study is concerned with the near-field region (i.e., the field region near the 
origin of the jet), whereas other studies, such as Gardiner & Prank (2000), are concerned 
with the far-field region. 

There is currently some debate over the origin of the knots in protostellar flows. It 
remains unclear as to whether the knots are due to the pulsing of the source of the jet or to 
instabilities in the beam. Our results show that for a given constant ejection velocity from 
the disk surface which produces a steady state jet, an episodic jet is produced when the 
given ejection velocity is pulsed. We showed that the temporal and spatial knotty structure 
of the jets depends on the period of the velocity perturbation at the base of the jet. 

Ouyed & Pudritz (1997b, 1999) rule out the possibility of a pinch MHD instabihty 
causing the knotty structure. This instability depends upon the strength of f3 and no such 
dependence was found by the authors. They found that only for /?j around unity, were knots 
produced. Birkinshaw (1990) noted that the Kelvin- Helmholtz (KH) instabihty does not 
exist if the beam speed is less than the Alfven speed. Hardee & Rosen (1999) showed that 
KH growth rates increase as the magnetosonic Mach number decreases, provided the jet is 
super- Alfvenic and the magnetized jet is nearly stabilized to the KH instability when the 
jet is sub- Alfvenic. Hardee & Rosen (2002) found that the jets can be stalilized to the KH 
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helical and higher order asymmetric normal modes provided the velocity shear, = Uj — Ue 
between the jet and the external medium is less than a "surface" Alfven speed, vas = 
[{Pj + Pe) {Bj + B'^) I (47rpjp£;)]^^^, where Pj{pe) is the jet (external medium) density and 
Bj{Be) is the jet (external medium) magnetic field. Ouyed et al. (2003) showed that the 
jet beyound the Alfven surface becomes unstable to nonaxisymmetric KH instabilities, in 
three-dimensional simulations of jets from Keplerian accretion disks. 

We also made a number of simulations using values of (5i different from unity. Figure 
11 shows the results of one of these simulations, a periodic perturbation at the disk surface, 
using a period of T = 20. We plotted the isocontours of the density. The results for j3i — 0.1 
are shown at the top of the figure. For this case, we see a smoothing of the structures along 
the jet. The isocontours of the density oscillate shghtly near the jet axis. At the botton of 
this figure, the results for /3j = 10 are shown. Here we observe the formation of structures 
from the jet axis up to a distance of about z — SSr^. The region near the jet axis shows 
irregular structures and is less well defined (compare with Figure 5, where (5i — l). 

Our results are consistent with those of Hardee & Rosen (1999, 2002) and Ouyed et 
al. (2003). This is shown in the sequence of simulations in Figures 5 and 11. When /3j is 
small (strong magnetic fields, and the flow is sub-Alfvenic), there is a not much growth of 
the instability. The super- Alfvenic case in which f3i is close to unity (the flow speed is close 
to the Alfven speed) is more unstable. 

EDITOR: PLACE FIGURE 11 HERE. 

We simulated a region very near the central object (20ri x 80^). It is interesting to 
evaluate the spatial dimensions involved in our simulations. For a protostar of mass O.SM©, 
taking = 0.04 AU for the inner radius of the disk, we have SOr^ ~ 10~^ pc. For a black hole 
of 1O^M0 (the Schwartzschild radius is ~ 1 AU), assuming = 20.6 AU, we have SOr^ ~ 10~^ 
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pc. A distance of separation of ~ llr^ of the structures corresponds, for the protostar case, to 
~ 10~^ pc and for the black hole case, to ~ 10~^ pc. These distances are much smaller than 
those between the Herbig-Haro objects observed in protostellar jets (about 10~^ pc). They 
are also much smaller than the distances observed between the emission knots, frequently 
found in AGN jets. In contrast to other simulations (e.g.. Prank et al.(2000)), we obtain 
structures very near the formation region of the jets, involving distances not yet resolved by 
current instruments. 

We can not say, at the present moment, whether or not the origin of the knotty structures 
seen in many jets at large distances from the source are due to the structures observed in the 
present simulations for regions close to the source. Further studies are needed to investigate 
this question. 

The instabilities seen by Ouyed & Pudritz (1997b, 1999), producing the knotty struc- 
tures, occur at the fast magnetosonic surface and propagate downstream. These instabilities, 
as explained by the authors, are due to MHD shocks. Some conditions in the knot generating 
region (in particular, for high B^/B^) seem to point towards a scenario in which the knot 
spacings could possibly be identified with the resonant wavelength pinch mode. As noted 
above, Ouyed & Pudritz (1997b, 1999) rule out the possibility of a pinch MHD instability 
causing the knotty structures since it has a strong /? dependence, which was not found by 
the authors. The authors only observed knotty structures when /3 was around unity. 

We believe that the knotty structures that we obtained in our simulations are due to the 
association of the pulsing mass load at the base of the jet and the instability due to MHD 
shocks, discussed in the studies of Ouyed & Pudritz (1997b, 1999). Our results indicate that 
a particular frequency of the pulsed velocity can excite the instability and produce persistent 
knotty structures along the jet. 
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Fig. 1. — Setup of the numerical simulations. Ghost zones {z < 0) define the surface of 
the Keplerian disk, while active zones (which are evolved in time) define the corona. The 
poloidal magnetic field continues smoothly into the disk surface (Fig. 2 of OP97a). 

Fig. 2. — Initial conditions in the atmosphere: density distribution in units of pi (eq. [8]) at 
the top and the magnetic potential (eq. [10]) at the bottom. Forty logarithmically spaced 
isocontour lines are shown for the density; 40 linearly spaced contours for the magnetic 
potential. The absolute nominal value for the steps between the isocontours of the logarithm 
of the density is 0.098. The velocity and toroidal magnetic field are initially zero. 

Fig. 3. — Simulation A (without perturbation): (from top to bottom) Density, poloidal 
velocity, poloidal magnetic field, toroidal magnetic field and toroidal velocity for r = 400. 
Forty logarithmically spaced isocontour lines are shown for the density, and 40 linearly spaced 
contours for all the other quantities. The labels of the toroidal magnetic field and toroidal 
velocity denote 100 times their nominal values, in units of and vkz, respectively. The 
absolute nominal values for the steps between the isocontours of the logarithm of the density, 
toroidal magnetic field and toroidal velocity are 0.1, 0.034Sj and O.OllSvKi: respectively. The 
maximum value of the poloidal velocity is equal to 0.64^^^:^. We note the formation of a steady 
state jet. 

Fig. 4. — Simulation B (T = 10): (from top to bottom) Density, poloidal velocity, poloidal 
magnetic field, toroidal magnetic field and toroidal velocity for r = 400. Forty logarithmically 
spaced isocontour lines are shown for the density, and 40 linearly spaced contours for all the 
other quantities. The labels of the toroidal magnetic field and toroidal velocity denote 
100 times their nominal values. The absolute nominal values for the steps between the 
isocontours of the logarithm of the density, toroidal magnetic field and toroidal velocity are 
0.08, 0.0368Sj and O.OlOSvKi, respectively. The maximum value of the poloidal velocity is 
equal to O.QivKi- We note the formation of regularly spaced structures along the jet axis. 
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Fig. 5. — Simulation C {T — 20): (from top to bottom) Density, poloidal velocity, poloidal 
magnetic field, toroidal magnetic field and toroidal velocity for r = 400, using a better 
resolution (250, 80) uniform zones. Forty logarithmically spaced isocontour lines are shown 
for the density, and 40 linearly spaced contours for all the other quantities. The labels 
of the toroidal magnetic field and toroidal velocity denote 100 times their nominal values. 
The absolute nominal values for the steps between the isocontours of the logarithm of the 
density, toroidal magnetic field and toroidal velocity are 0.094, 0.0346i?j and 0.0154vKi, 
respectively. The maximum value of the poloidal velocity is equal to 0.61vKi- We note, 
again, the formation of regularly spaced structures along the jet axis. 

Fig. 6. — Simulation D (T — 40): (from top to bottom) Density, poloidal velocity, poloidal 
magnetic field, toroidal magnetic field and toroidal velocity for r = 400. Forty logarithmically 
spaced isocontour lines are shown for the density, and 40 linearly spaced contours for all the 
other quantities. The labels of the toroidal magnetic field and toroidal velocity denote 
100 times their nominal values. The absolute nominal values for the steps between the 
isocontours of the logarithm of the density, toroidal magnetic field and toroidal velocity are 
0.094, 0.0342Sj and O.OldvKi, respectively. The maximum value of the poloidal velocity is 
equal to 0.67vKi- We note fragmentation of regularly spaced structures along the jet axis. 

Fig. 7. — Simulation E {T — 80): (from top to bottom) Density, poloidal velocity, poloidal 
magnetic field, toroidal magnetic field and toroidal velocity for r = 400. Forty logarithmically 
spaced isocontour lines are shown for the density, and 40 linearly spaced contours for all the 
other quantities. The labels of the toroidal magnetic field and toroidal velocity denote 
100 times their nominal values. The absolute nominal values for the steps between the 
isocontours of the logarithm of the density, toroidal magnetic field and toroidal velocity are 
0.096, O.OSMBi and 0.0152vKi, respectively. The maximum value of the poloidal velocity is 
equal to O.dAvxi- We note more pronounced fragmentation of the regularly spaced structures 
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along the jet axis. 

Fig. 8. — Simulations A to E: Density (solid line), axial velocity (dotted line), toroidal 
velocity (dashed hne) and toroidal magnetic field (dashed-dotted line) along the 2;-axis at 
r = Srj. The curves are normalized by their maximum values. 

Fig. 9. — Simulations A to E: Axial and radial mass fiuxes (eq.[18]). Here we show the axial 
(solid line) and radial (dotted line) mass fluxes at r = 400. We calculate the axial mass flux, 
at each axial distance using the transverse area A — irr^ with r = 20ri. For the radial 
mass flux, we use the cylindrical surface area A — 27rrz with r = 20ri, at each axial distance 
z. The fluxes are given in units of rrii = 2'nr'jpiVKi- The scale on the vertial axis represents 
a tenth of the nominal value. 

Fig. 10. — Simulations A to E: The sonic Mach numbers (dotted lines), Alfvenic Mach 
numbers (dashed-dotted lines) and magnetosonic Mach numbers (solid lines) along the z- 
axis for all simulations at r = 400 in a cut at r = 5rj. For simulatations C to E, dashed-dotted 
lines denote 10~^ times the nominal values of Alfvenic Mach numbers. 

Fig. 11. — Logarithmically spaced isocontour lines are shown for the density. On the top 
we have the result Pi = 0.1. We obtained, in this case, smoothing of the structures along 
the jet. The isocontours of the density shghtly oscillate near the jet axis. On the botton of 
the figure we used /3j = 10. Here we observe the formation of structures displaced from the 

jet axis up to the distance of about z = SSr^. The region near the jet axis show irregular 
structures and is less well defined. 



